{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "178b8bcd",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62207426",
   "metadata": {},
   "source": [
    "# Домашнее задание №12\n",
    "\n",
    "# 1 Задание №1\n",
    "\n",
    "На любом языке программирования напишите программу, вычисляющую по заданной\n",
    "таблице сопряжённости признаков\n",
    "\n",
    "- статистику хи-квадрат,\n",
    "- P-значение для критерия хи-квадрат Фишера-Пирсона,\n",
    "- нормированные коэффициенты связи, основанные на статистике хи-квадрат,\n",
    "- меры прогнозы Гутмана,\n",
    "- коэффициенты контингенции и ассоциации (если таблица имеет размер 2 × 2).\n",
    "\n",
    "## 1.1 Подготовка функций для решения задания №1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "467f2612",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_p_value(stat_value, side='both', stat_type='norm', df=1):\n",
    "    \"\"\" Расчет p-value для значения статистики\n",
    "    \n",
    "    Параметры\n",
    "    ---------\n",
    "    stat_value : float\n",
    "      Значение статистики\n",
    "    side : str\n",
    "      Сторона теста для вычисления p-value\n",
    "    stat_type : str\n",
    "      Тип распределения значения статистики.\n",
    "      Поддерживаемые значения: norm, ch2 и t\n",
    "    df : int\n",
    "      Количество степеней свободы. Только для chi2 и t\n",
    "    Результат\n",
    "    ---------\n",
    "    p_val : float\n",
    "    \"\"\"\n",
    "    \n",
    "    side = str(side).lower().strip()\n",
    "    assert side in ('both', 'left', 'right')\n",
    "    stat_type = str(stat_type).lower().strip()\n",
    "    assert stat_type in ('norm', 'chi2', 't')\n",
    "    \n",
    "    if stat_type == 'norm':\n",
    "        stat_quantile = stats.norm.cdf(stat_value)\n",
    "    elif stat_type == 'chi2':\n",
    "        stat_quantile = stats.chi2.cdf(stat_value, df)\n",
    "    else:\n",
    "        stat_quantile = stats.t.cdf(stat_value, df)\n",
    "    \n",
    "    if side == 'both':\n",
    "        p_val = np.min([2*stat_quantile, 2 - 2*stat_quantile])\n",
    "    elif side == 'left':\n",
    "        p_val = stat_quantile\n",
    "    else:\n",
    "        p_val = 1 - stat_quantile\n",
    "        \n",
    "    return p_val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "546d0581",
   "metadata": {},
   "outputs": [],
   "source": [
    "def chi2_test(array, return_stats_norm=False):\n",
    "    \"\"\"Расчет статистики Хи-квадрат\n",
    "    \n",
    "    Параметры\n",
    "    ---------\n",
    "    array : np.ndarray, pd.DataFrame\n",
    "      Массив данных\n",
    "    return_stats_norm : bool\n",
    "      Если True, будут расчитаны дополнительные нормированные \n",
    "      коэффициенты статистики\n",
    "      \n",
    "    Результат\n",
    "    ---------\n",
    "    Если return_stats_norm == False:\n",
    "    chi2, p_value\n",
    "    \n",
    "    иначе:\n",
    "    chi2, p_value, chi2_norm\n",
    "    \n",
    "    chi2 : float\n",
    "      Статистика\n",
    "    p_value : float\n",
    "      p-p_value для статистики\n",
    "    chi2_norm : dict\n",
    "      Словарь нормированных коэффициентов статистики:\n",
    "      {\n",
    "        # Коэффициент Пирсона\n",
    "        'pearson': float,\n",
    "        # Коэффициент Чупрова\n",
    "        'chuprov': float,\n",
    "        # Коэффициент Крамера\n",
    "        'kramer': float,\n",
    "        # Коэффициент среднеквадратичной сопряженности\n",
    "        'sks': float,\n",
    "    }\n",
    "    \"\"\"\n",
    "    \n",
    "    if isinstance(array, pd.DataFrame):\n",
    "        array = array.to_numpy()\n",
    "        \n",
    "    if not isinstance(array, np.ndarray):\n",
    "        array = np.array(array)\n",
    "      \n",
    "    \n",
    "    if np.min(array) == 0:\n",
    "        print('ВНИМАНИЕ: Найдены нулевые частоты в таблице контенгенций, тест может работать некорректно')\n",
    "    elif np.min(array) <= 5:\n",
    "        print('ВНИМАНИЕ: Найдены частоты менее 5 в таблице контенгенций, тест может работать некорректно')\n",
    "        \n",
    "    n = np.sum(array)\n",
    "    # Сумма каждой строки отдельно\n",
    "    n_i = np.sum(array, axis=1)\n",
    "    # Сумма каждого столбца отдельно\n",
    "    n_j = np.sum(array, axis=0)\n",
    "    \n",
    "    # Кол-во строк\n",
    "    k = array.shape[0]\n",
    "    # Кол-во столбцов\n",
    "    m = array.shape[1]\n",
    "    \n",
    "    assert np.sum(n_i) == n\n",
    "    assert np.sum(n_j) == n\n",
    "    \n",
    "    chi2 = []\n",
    "    \n",
    "    for i in range(array.shape[0]):\n",
    "        temp = 0\n",
    "        for j in range(array.shape[1]):\n",
    "            a =  n_i[i] * n_j[j] / n\n",
    "            temp += np.square(array[i,j] - a ) / a  \n",
    "        chi2.append(temp)\n",
    "        \n",
    "    chi2 = np.sum(chi2)\n",
    "    \n",
    "    # Тест правосторонний\n",
    "    p_value = get_p_value(chi2, side='right', stat_type='chi2', df=(k-1) * (m-1))\n",
    "    \n",
    "    if not return_stats_norm:\n",
    "        return chi2, p_value\n",
    "    \n",
    "    chi2_norm = {\n",
    "        'pearson': np.sqrt(chi2 / (n + chi2)),\n",
    "        'chuprov': np.sqrt( chi2 / (n * np.sqrt( (k-1) * (m-1) )) ),\n",
    "        'kramer': np.sqrt( chi2 / (n * np.min( [(k-1), (m-1)] )) ),\n",
    "        'sks': np.sqrt( chi2 / n ),\n",
    "    }\n",
    "    \n",
    "    assert all([0<=stat<=1 for stat in chi2_norm.values()]), 'Найдены нормированные статистики вне диапазона от 0 до 1'\n",
    "    \n",
    "    return chi2, p_value, chi2_norm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "0a7b2e5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def gutman_test(array):\n",
    "    \"\"\"Мера прогноза Гутмана\n",
    "    \n",
    "    Параметры\n",
    "    ---------\n",
    "    array : np.ndarray, pd.DataFrame\n",
    "      Массив данных\n",
    "      \n",
    "    Результат\n",
    "    ---------\n",
    "    lambda : float\n",
    "    \"\"\"\n",
    "    \n",
    "    if isinstance(array, pd.DataFrame):\n",
    "        array = array.to_numpy()\n",
    "        \n",
    "    if not isinstance(array, np.ndarray):\n",
    "        array = np.array(array)\n",
    "     \n",
    "    lammdas = []\n",
    "    \n",
    "    n = np.sum(array)\n",
    "    \n",
    "    # Ищем вероятности ошибок прогнозов для обычной таблицы и транспонированной.\n",
    "    # Результаты будем усреднять.\n",
    "    for array_variant in (array, array.T):\n",
    "    \n",
    "        n_j = np.sum(array, axis=0)\n",
    "        \n",
    "        p_b1 = 1 - np.max(n_j) / n\n",
    "        p_b2 = 1 - np.sum(np.max(array, axis=1) / n)\n",
    "        \n",
    "        lambda_variant = (p_b1 - p_b2) / p_b1\n",
    "        lammdas.append(lambda_variant)\n",
    "        \n",
    "    return np.mean(lammdas)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "8c303a29",
   "metadata": {},
   "outputs": [],
   "source": [
    "def contingent_coef(array):\n",
    "    \"\"\" Коэффициент контингенции\n",
    "    \n",
    "    Параметры\n",
    "    ---------\n",
    "    array : np.ndarray, pd.DataFrame\n",
    "      Массив данных\n",
    "      \n",
    "    Результат\n",
    "    ---------\n",
    "    contingent_coef : float\n",
    "    \"\"\"\n",
    "    \n",
    "    if isinstance(array, pd.DataFrame):\n",
    "        array = array.to_numpy()\n",
    "        \n",
    "    if not isinstance(array, np.ndarray):\n",
    "        array = np.array(array)\n",
    "        \n",
    "    assert array.shape[0] == 2 and array.shape[1] == 2\n",
    "    \n",
    "    a, b, c, d = array.flatten()\n",
    "    \n",
    "    return (a*d - b*c) / np.sqrt( (a+b) * (c+d) * (a+c) * (b+d) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a09f66c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def association_coef(array):\n",
    "    \"\"\" Коэффициент ассоциации\n",
    "    \n",
    "    Параметры\n",
    "    ---------\n",
    "    array : np.ndarray, pd.DataFrame\n",
    "      Массив данных\n",
    "      \n",
    "    Результат\n",
    "    ---------\n",
    "    contingent_coef : float\n",
    "    \"\"\"\n",
    "    \n",
    "    if isinstance(array, pd.DataFrame):\n",
    "        array = array.to_numpy()\n",
    "        \n",
    "    if not isinstance(array, np.ndarray):\n",
    "        array = np.array(array)\n",
    "        \n",
    "    assert array.shape[0] == 2 and array.shape[1] == 2\n",
    "    \n",
    "    a, b, c, d = array.flatten()\n",
    "    \n",
    "    return (a*b - b*c) / (a*d + b*c)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07122d63",
   "metadata": {},
   "source": [
    "\n",
    "## 1.2 Решение задач в задании №1\n",
    "### 1.2.1 Задача №1.1\n",
    "\n",
    "По данным исследования «Мониторинг социальных и экономических перемен в России» составлена следующая таблица сопряженности:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "927c2cc6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Благоприятная</th>\n",
       "      <th>Напряженная</th>\n",
       "      <th>Взрывоопасная критическая</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Хорошее</th>\n",
       "      <td>12</td>\n",
       "      <td>48</td>\n",
       "      <td>47</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Среднее</th>\n",
       "      <td>20</td>\n",
       "      <td>478</td>\n",
       "      <td>6666</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Плохое</th>\n",
       "      <td>11</td>\n",
       "      <td>160</td>\n",
       "      <td>701</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         Благоприятная  Напряженная  Взрывоопасная критическая\n",
       "Хорошее             12           48                         47\n",
       "Среднее             20          478                       6666\n",
       "Плохое              11          160                        701"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.DataFrame({\n",
    "    'Благоприятная': [12, 20, 11],\n",
    "    'Напряженная': [48, 478, 160],\n",
    "    'Взрывоопасная критическая': [47, 6666, 701]\n",
    "}, index = ['Хорошее', 'Среднее', 'Плохое'])\n",
    "df "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ffa2eb7",
   "metadata": {},
   "source": [
    "**Вопрос А:** «Как бы Вы оценили в настоящее время материальное положение\n",
    "вашей семьи?» (Хорошее, среднее, плохое).\n",
    "\n",
    "**Вопрос В:** «Как бы Вы оценили в целом политическую обстановку в России?»\n",
    "(Благополучная, напряжённая, взрывоопасная критическая).\n",
    "\n",
    "Имеется ли зависимость между признаками А и В? Прокомментируйте характер\n",
    "связи между А и В с помощью различных коэффициентов.\n",
    "\n",
    "### Решение\n",
    "\n",
    "Воспользуемся функцией chi2_test для получения всех необходимых данных:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "5aecdfd1",
   "metadata": {},
   "outputs": [],
   "source": [
    "chi2, p_value, chi2_norm = chi2_test(df, return_stats_norm=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "6d7b3f65",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "====================================\n",
      "chi2: 585.4980121554853\n",
      "p_value: 0.0\n",
      "pearson: 0.258995938581778\n",
      "chuprov: 0.18960749353865317\n",
      "kramer: 0.18960749353865317\n",
      "sks: 0.2681454888899323\n",
      "guttman: 0.0013717421124842277\n",
      "====================================\n",
      "Между признаками имеется зависимость на уровне значимости 0.05\n"
     ]
    }
   ],
   "source": [
    "alpha = 0.05\n",
    "\n",
    "print(\"====================================\")\n",
    "print('chi2:', chi2)\n",
    "print('p_value:', p_value)\n",
    "for chi2_norm_coef_name, chi2_norm_coef_val in chi2_norm.items():\n",
    "    print(f\"{chi2_norm_coef_name}: {chi2_norm_coef_val}\")\n",
    "print('guttman:', gutman_test(df))\n",
    "\n",
    "print(\"====================================\")\n",
    "if p_value < alpha:\n",
    "    print(f\"Между признаками имеется зависимость на уровне значимости {alpha}\")\n",
    "else:\n",
    "    print(f\"Между признаками отсутствует зависимость на уровне значимости {alpha}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "856746f2",
   "metadata": {},
   "source": [
    "Тест определил наличие зависимости. P-value равно 0, поэтому с точки зрения теста вероятность зависимости очень высокая даже при меньшем пороге уровне значимости. При этом нормированные коэффициенты говорят о наличии очень слабой зависимости, так как они не достигают даже уровня 0.5\n",
    "\n",
    "### 1.2.2 Задача №1.2\n",
    "\n",
    "Изучается зависимость выбираемой абитуриентами специальности от пола. Было опрошено 480 человек и получены следующие результаты:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "24204f1e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>М</th>\n",
       "      <th>Ж</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Естественные науки</th>\n",
       "      <td>168</td>\n",
       "      <td>85</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Гуманитарные науки</th>\n",
       "      <td>85</td>\n",
       "      <td>135</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                      М    Ж\n",
       "Естественные науки  168   85\n",
       "Гуманитарные науки   85  135"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.DataFrame({\n",
    "    'М': [168, 85],\n",
    "    'Ж': [85, 135],\n",
    "}, index = ['Естественные науки', 'Гуманитарные науки'])\n",
    "df "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90f357f8",
   "metadata": {},
   "source": [
    "Можно ли на уровне доверия 0,95 считать, что между полом и выбором специальности есть связь? Оцените силу связи с помощью различных коэффициентов."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "8de1c82e",
   "metadata": {},
   "outputs": [],
   "source": [
    "chi2, p_value, chi2_norm = chi2_test(df, return_stats_norm=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "a61d1cce",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "====================================\n",
      "chi2: 36.468067967004636\n",
      "p_value: 1.5518788476498457e-09\n",
      "pearson: 0.26754564947027903\n",
      "chuprov: 0.2776679841897233\n",
      "kramer: 0.2776679841897233\n",
      "sks: 0.2776679841897233\n",
      "guttman: 0.22727272727272743\n",
      "contingent_coef: 0.2776679841897233\n",
      "association_coef: 0.23591372680153821\n",
      "====================================\n",
      "Между признаками имеется зависимость на уровне значимости 0.05\n"
     ]
    }
   ],
   "source": [
    "alpha = 0.05\n",
    "\n",
    "print(\"====================================\")\n",
    "print('chi2:', chi2)\n",
    "print('p_value:', p_value)\n",
    "for chi2_norm_coef_name, chi2_norm_coef_val in chi2_norm.items():\n",
    "    print(f\"{chi2_norm_coef_name}: {chi2_norm_coef_val}\")\n",
    "print('guttman:', gutman_test(df))\n",
    "print('contingent_coef:', contingent_coef(df))\n",
    "print('association_coef:', association_coef(df))\n",
    "\n",
    "\n",
    "print(\"====================================\")\n",
    "if p_value < alpha:\n",
    "    print(f\"Между признаками имеется зависимость на уровне значимости {alpha}\")\n",
    "else:\n",
    "    print(f\"Между признаками отсутствует зависимость на уровне значимости {alpha}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a4c66b4",
   "metadata": {},
   "source": [
    "Аналогично тест определил наличие зависимости. И точно также нормированные коэффициенты говорят  о наличии очень слабой зависимости, так как они не достигают даже уровня 0.5\n",
    "\n",
    "# 2 Задание №2\n",
    "\n",
    "На любом языке программирования напишите программу, которая по результатам наблюдений за двумя количественными признаками:\n",
    "- вычисляет выборочный коэффициент корреляции и P-значение для критерия проверки гипотезы о нулевом значении коэффициента корреляции;\n",
    "- вычисляет коэффициент Спирмена и P-значение для критерия Спирмена;\n",
    "- вычисляет коэффициент Крамера и P-значение для критерия Крамера;\n",
    "- разбивая диапазон количественных признаков на несколько интервалов строит по наблюдениям таблицу сопряжённости, вычисляет нормированные коэффициенты связи (Пирсона, Крамера и др.), меры прогноза Гутмана и P-значение для критерия хи-квадрат Фишера Пирсона.\n",
    "\n",
    "## 2.1 Подготовка функций для решения задания №2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "82985cbf",
   "metadata": {},
   "outputs": [],
   "source": [
    "def pearson_test(array):\n",
    "    \"\"\" Тест Пирсона (стандартный корреляционный)\n",
    "    \n",
    "    Можно применять только для данных с двумя наборами признаков.\n",
    "    \n",
    "    Параметры\n",
    "    ---------\n",
    "    array : np.ndarray, pd.DataFrame\n",
    "      Массив данных\n",
    "      \n",
    "    Результат\n",
    "    ---------\n",
    "    stat, p_value, r_xy : tuple\n",
    "    \n",
    "    stat : float\n",
    "      Значение статистики\n",
    "    p_value : float\n",
    "      p-value для статистики\n",
    "    r_xy : float\n",
    "      Коэффициент корреляции\n",
    "    \"\"\"\n",
    "    \n",
    "    if isinstance(array, pd.DataFrame):\n",
    "        array = array.to_numpy()\n",
    "        \n",
    "    if not isinstance(array, np.ndarray):\n",
    "        array = np.array(array)\n",
    "        \n",
    "    assert array.shape[1] == 2\n",
    "    \n",
    "    n = array.shape[0] \n",
    "    \n",
    "    means = array.mean(axis=0)\n",
    "    array_centered = array - means\n",
    "    \n",
    "    r_xy = np.sum(array_centered[:, 0] * array_centered[:, 1]) / np.sqrt( np.sum(array_centered[:, 0]**2) \n",
    "                                                                          * np.sum(array_centered[:, 1]**2) )\n",
    "    assert -1 <= r_xy <= 1\n",
    "    \n",
    "    stat = r_xy * np.sqrt(n - 2) / np.sqrt(1 - r_xy**2)\n",
    "    p_value = get_p_value(stat, side='both', stat_type='t', df=n-2)\n",
    "    \n",
    "    return stat, p_value, r_xy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "eafb4e2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def spearman_test(array):\n",
    "    \"\"\" Тест Спирмана\n",
    "    \n",
    "    Можно применять только для данных с двумя наборами признаков.\n",
    "    \n",
    "    Параметры\n",
    "    ---------\n",
    "    array : np.ndarray, pd.DataFrame\n",
    "      Массив данных\n",
    "      \n",
    "    Результат\n",
    "    ---------\n",
    "    stat, p_value, p_s : tuple\n",
    "    \n",
    "    stat : float\n",
    "      Значение статистики\n",
    "    p_value : float\n",
    "      p-value для статистики\n",
    "    p_s : float\n",
    "      Коэффициент корреляции Спирмана\n",
    "    \"\"\"\n",
    "    \n",
    "    if isinstance(array, pd.DataFrame):\n",
    "        array = array.to_numpy()\n",
    "        \n",
    "    if not isinstance(array, np.ndarray):\n",
    "        array = np.array(array)\n",
    "        \n",
    "    assert array.shape[1] == 2\n",
    "    \n",
    "    n = array.shape[0]\n",
    "    \n",
    "    # Ранжируем каждый признак отдельно без сортировки\n",
    "    rank_data_x = stats.rankdata(array[:, 0])\n",
    "    rank_data_y = stats.rankdata(array[:, 1])\n",
    "    \n",
    "    # Находим сумму квадратов разниц рангов\n",
    "    s = np.sum( ( rank_data_x - rank_data_y)**2 )\n",
    "    \n",
    "    # Коэффициент корреляции Спирмана\n",
    "    p_s = 1 - 6*s / (n**3 - n)\n",
    "    \n",
    "    assert -1 <= p_s <= 1\n",
    "    \n",
    "    stat = np.sqrt(n - 1) * p_s\n",
    "    p_value = get_p_value(stat, side='both', stat_type='norm')\n",
    "    \n",
    "    return stat, p_value, p_s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "012c77fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "def kendall_test(array):\n",
    "    \"\"\" Тест Кендалла\n",
    "    \n",
    "    Можно применять только для данных с двумя наборами признаков.\n",
    "    \n",
    "    Параметры\n",
    "    ---------\n",
    "    array : np.ndarray, pd.DataFrame\n",
    "      Массив данных\n",
    "      \n",
    "    Результат\n",
    "    ---------\n",
    "    stat, p_value, p_s : tuple\n",
    "    \n",
    "    stat : float\n",
    "      Значение статистики\n",
    "    p_value : float\n",
    "      p-value для статистики\n",
    "    t : float\n",
    "      Коэффициент корреляции Кендалла\n",
    "    \"\"\"\n",
    "    \n",
    "    \n",
    "    if isinstance(array, pd.DataFrame):\n",
    "        array = array.to_numpy()\n",
    "        \n",
    "    if not isinstance(array, np.ndarray):\n",
    "        array = np.array(array)\n",
    "        \n",
    "    assert array.shape[1] == 2\n",
    "    \n",
    "    n = array.shape[0]\n",
    "    \n",
    "    # Ранжируем каждый признак отдельно без сортировки\n",
    "    rank_data_x = stats.rankdata(array[:, 0])\n",
    "    rank_data_y = stats.rankdata(array[:, 1])\n",
    "    \n",
    "    # Получаем список пар рангов признаков [(rank_x1, rank_y1), (rank_x2, rank_y2), ...]\n",
    "    rank_pairs = list(zip(rank_data_x, rank_data_y))\n",
    "    \n",
    "    # Сортируем список пар по признаку x от меньшего к большему\n",
    "    sorted_rank_data = list(sorted(rank_pairs))\n",
    "    \n",
    "    # В итоге получаем отсортированные ранги признака Y по признаку X\n",
    "    sorted_rank_data = np.array(list(zip(*sorted_rank_data)))\n",
    "    # Для нашей формулы требуются только ранги признака Y\n",
    "    sorted_rank_data_y = sorted_rank_data[1]\n",
    "    \n",
    "    # Кол-во элементов, у которых ранг выше текущей позиции\n",
    "    y_up = 0\n",
    "    # Кол-во элементов, у которых ранг ниже текущей позиции\n",
    "    y_down = 0\n",
    "    \n",
    "    for idx in range(n-1):\n",
    "        y_up += len(np.where(sorted_rank_data_y[idx+1:] > sorted_rank_data_y[idx])[0])\n",
    "        y_down += len(np.where(sorted_rank_data_y[idx+1:] < sorted_rank_data_y[idx])[0])  \n",
    "        \n",
    "    s = y_up - y_down\n",
    "    t = 2 * s / (n * (n - 1))\n",
    "    \n",
    "    assert -1 <= t <= 1\n",
    "    \n",
    "    stat = 3 * np.sqrt(n) * t / 2\n",
    "    p_value = get_p_value(stat, side='both', stat_type='norm')\n",
    "    \n",
    "    return stat, p_value, t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "da0a9eb1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_contingency_table(array, bins=2):\n",
    "    \"\"\" Формирование таблицы сопряженности\n",
    "    \n",
    "    Функция разбивает на интервалы каждый признак по отдельности.\n",
    "    Поддерживаются только массивы с двумя числовыми признаками.\n",
    "    Формируется матрица, у которой количество строк - количество бинов первого признака,\n",
    "    а количество столбцов - количество бинов второго признака.\n",
    "    Берется значение с двумя признаками и счётчик увеличивается в той ячейке,\n",
    "    в диапазон которых попадают оба значения.\n",
    "\n",
    "    Параметры\n",
    "    ---------\n",
    "    array : np.ndarray, pd.DataFrame\n",
    "      Массив данных\n",
    "    bins : int\n",
    "      Количество интервалов, на которые нужно разбить данные\n",
    "      \n",
    "    Результат\n",
    "    ---------\n",
    "    contingency_table : np.ndarray\n",
    "      Имеет размер B x C,\n",
    "      где С - кол-во столбцов в array и B - количество бинов\n",
    "    \"\"\"\n",
    "    \n",
    "    if isinstance(array, pd.DataFrame):\n",
    "        array = array.to_numpy()\n",
    "        \n",
    "    if not isinstance(array, np.ndarray):\n",
    "        array = np.array(array)\n",
    "        \n",
    "    assert array.shape[1] == 2\n",
    "    \n",
    "    # Здесь получаем интервалы для каждого отдельного признака: [интервалы_признак_1, интервалы_признак_2, .. ]\n",
    "    categories_ranges = [pd.cut(features, bins=bins).categories for features in array.T]\n",
    "    \n",
    "    # Делаем таблицу сопряженности с нулевыми значениями размером bins x bins\n",
    "    contingency_table = np.zeros(shape=[bins, bins])\n",
    "    \n",
    "    # Пробегаемся по каждой колонке и строке\n",
    "    for row_idx in range(array.shape[0]):\n",
    "        # Значение исходной таблицы\n",
    "        value1 = array[row_idx, 0]\n",
    "        value2 = array[row_idx, 1]\n",
    "        # Находим к какому номеру диапазона относится значение\n",
    "        category_1 = np.argmax([value1 in value_range for value_range in categories_ranges[0]])\n",
    "        # Находим к какому номеру диапазона относится значение\n",
    "        category_2 = np.argmax([value2 in value_range for value_range in categories_ranges[1]])\n",
    "        # Увеличиваем значение на 1 для нужного номера диапазона\n",
    "        contingency_table[category_1, category_2] += 1\n",
    "        \n",
    "    # Удаляем все строки, у которых только нулевые значения\n",
    "    contingency_table = contingency_table[contingency_table.any(axis=1)]\n",
    "    # Удаляем все столбцы, у которых только нулевые значения\n",
    "    contingency_table = contingency_table[:, contingency_table.any(axis=0)]\n",
    "             \n",
    "    return contingency_table"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b509dc3",
   "metadata": {},
   "source": [
    "## 2.2 Решение задач в задании №2\n",
    "### 2.2.1 Задача №2.1\n",
    "\n",
    "Группа пловцов из 15 человек принимает участие в местной спортивной олимпиаде. В программе состязаний два заплыва: 50 м вольным стилем и 100 м баттерфляем. \n",
    "\n",
    "Результаты соревнований следующие:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "b7bef727",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>50 м вольн.</th>\n",
       "      <th>100 м батт.</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>22.49</td>\n",
       "      <td>52.93</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>22.56</td>\n",
       "      <td>53.40</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>23.45</td>\n",
       "      <td>53.70</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>22.58</td>\n",
       "      <td>53.36</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>24.30</td>\n",
       "      <td>61.80</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>24.20</td>\n",
       "      <td>55.20</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>23.47</td>\n",
       "      <td>53.54</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>23.50</td>\n",
       "      <td>58.33</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>24.48</td>\n",
       "      <td>60.40</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>25.02</td>\n",
       "      <td>60.30</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>23.04</td>\n",
       "      <td>54.28</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>23.24</td>\n",
       "      <td>53.60</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>25.20</td>\n",
       "      <td>62.24</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>24.61</td>\n",
       "      <td>54.45</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>26.02</td>\n",
       "      <td>61.52</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    50 м вольн.  100 м батт.\n",
       "1         22.49        52.93\n",
       "2         22.56        53.40\n",
       "3         23.45        53.70\n",
       "4         22.58        53.36\n",
       "5         24.30        61.80\n",
       "6         24.20        55.20\n",
       "7         23.47        53.54\n",
       "8         23.50        58.33\n",
       "9         24.48        60.40\n",
       "10        25.02        60.30\n",
       "11        23.04        54.28\n",
       "12        23.24        53.60\n",
       "13        25.20        62.24\n",
       "14        24.61        54.45\n",
       "15        26.02        61.52"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = [22.49, 22.56, 23.45, 22.58, 24.3, 24.2, 23.47, 23.5, 24.48, 25.02, 23.04, 23.24, 25.2, 24.61, 26.02]\n",
    "b = [52.93, 53.4, 53.7, 53.36, 61.8, 55.2, 53.54, 58.33, 60.4, 60.3, 54.28, 53.6, 62.24, 54.45, 61.52]\n",
    "\n",
    "df = pd.DataFrame([a, b], columns=np.arange(1, len(a)+1), index=['50 м вольн.', '100 м батт.']).transpose()\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b3401c4",
   "metadata": {},
   "source": [
    "Определите степень связи между результатами в разных заплывах."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "c337eb83",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "====================================\n",
      "\n",
      "pearson_stat: 4.9431604784651\n",
      "pearson_p_value: 0.0002686773575546475\n",
      "pearson_corr: 0.8079169977411043 \n",
      "\n",
      "spearman_stat: 3.300676337618441\n",
      "spearman_p_value: 0.0009645208069288813\n",
      "spearman_corr: 0.8821428571428571 \n",
      "\n",
      "kendall_stat: 4.038968346759163\n",
      "kendall_p_value: 5.3686808705766254e-05\n",
      "kendall_corr: 0.6952380952380952 \n",
      "\n",
      "====================================\n",
      "\n",
      "Таблица сопряженности для 2 бинов:\n",
      "\n",
      "[[8. 1.]\n",
      " [1. 5.]] \n",
      "\n",
      "ВНИМАНИЕ: Найдены частоты менее 5 в таблице контенгенций, тест может работать некорректно\n",
      "chi2: 7.824074074074074\n",
      "p_value: 0.005155485148346517\n",
      "pearson: 0.5854905538443584\n",
      "chuprov: 0.7222222222222222\n",
      "kramer: 0.7222222222222222\n",
      "sks: 0.7222222222222222\n",
      "guttman: 7.484210060072233e-16\n",
      "\n",
      "Между признаками в таблице сопряженности имеется зависимость на уровне значимости 0.05\n"
     ]
    }
   ],
   "source": [
    "alpha = 0.05\n",
    "bins = 2\n",
    "\n",
    "print(\"====================================\\n\")\n",
    "for test_name, test_func in zip(['pearson', 'spearman', 'kendall'], [pearson_test, spearman_test, kendall_test]):\n",
    "    stat, p_value, corr = test_func(df)\n",
    "    print(f'{test_name}_stat:', stat)\n",
    "    print(f'{test_name}_p_value:', p_value)\n",
    "    print(f'{test_name}_corr:', corr, '\\n')\n",
    "\n",
    "contingency_table = make_contingency_table(df, bins=bins)\n",
    "\n",
    "print(\"====================================\\n\")\n",
    "print(f'Таблица сопряженности для {bins} бинов:\\n')\n",
    "print(contingency_table, '\\n')\n",
    "\n",
    "chi2, p_value, chi2_norm = chi2_test(contingency_table, return_stats_norm=True)\n",
    "\n",
    "print('chi2:', chi2)\n",
    "print('p_value:', p_value)\n",
    "for chi2_norm_coef_name, chi2_norm_coef_val in chi2_norm.items():\n",
    "    print(f\"{chi2_norm_coef_name}: {chi2_norm_coef_val}\")\n",
    "print('guttman:', gutman_test(df))\n",
    "\n",
    "if p_value < alpha:\n",
    "    print(f\"\\nМежду признаками в таблице сопряженности имеется зависимость на уровне значимости {alpha}\")\n",
    "else:\n",
    "    print(f\"\\nМежду признаками в таблице сопряженности отсутствует зависимость на уровне значимости {alpha}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c905bed",
   "metadata": {},
   "source": [
    "### 2.2.2 Задача №2.2\n",
    "\n",
    "В таблице приведены данные об урожайности пшеницы и картофеля на соседних полях. Что можно сказать о зависимости показателей урожайности картофеля и пшеницы?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "44108d89",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Пшеница</th>\n",
       "      <th>Картофель</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1926</th>\n",
       "      <td>20.1</td>\n",
       "      <td>7.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1927</th>\n",
       "      <td>23.6</td>\n",
       "      <td>7.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1928</th>\n",
       "      <td>26.3</td>\n",
       "      <td>7.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1929</th>\n",
       "      <td>19.9</td>\n",
       "      <td>6.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1930</th>\n",
       "      <td>16.7</td>\n",
       "      <td>6.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1931</th>\n",
       "      <td>23.2</td>\n",
       "      <td>7.3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1932</th>\n",
       "      <td>31.4</td>\n",
       "      <td>9.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1933</th>\n",
       "      <td>33.5</td>\n",
       "      <td>9.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1934</th>\n",
       "      <td>28.2</td>\n",
       "      <td>8.8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1935</th>\n",
       "      <td>35.3</td>\n",
       "      <td>10.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1936</th>\n",
       "      <td>29.3</td>\n",
       "      <td>8.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      Пшеница  Картофель\n",
       "1926     20.1        7.2\n",
       "1927     23.6        7.1\n",
       "1928     26.3        7.4\n",
       "1929     19.9        6.1\n",
       "1930     16.7        6.0\n",
       "1931     23.2        7.3\n",
       "1932     31.4        9.4\n",
       "1933     33.5        9.2\n",
       "1934     28.2        8.8\n",
       "1935     35.3       10.4\n",
       "1936     29.3        8.0"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = [20.1, 23.6, 26.3, 19.9, 16.7, 23.2, 31.4, 33.5, 28.2, 35.3, 29.3]\n",
    "b = [7.2, 7.1, 7.4, 6.1, 6.0, 7.3, 9.4, 9.2, 8.8, 10.4, 8.0]\n",
    "\n",
    "df = pd.DataFrame([a, b], columns=np.arange(1926, 1937), index=['Пшеница', 'Картофель']).transpose()\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "8d169496",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "====================================\n",
      "\n",
      "pearson_stat: 8.99104946703741\n",
      "pearson_p_value: 8.607825483286646e-06\n",
      "pearson_corr: 0.9485888243336862 \n",
      "\n",
      "spearman_stat: 3.0185377665243625\n",
      "spearman_p_value: 0.002539977411205818\n",
      "spearman_corr: 0.9545454545454546 \n",
      "\n",
      "kendall_stat: 4.251309958546467\n",
      "kendall_p_value: 2.1252380079328503e-05\n",
      "kendall_corr: 0.8545454545454545 \n",
      "\n",
      "====================================\n",
      "\n",
      "Таблица сопряженности для 2 бинов:\n",
      "\n",
      "[[5. 0.]\n",
      " [2. 4.]] \n",
      "\n",
      "ВНИМАНИЕ: Найдены нулевые частоты в таблице контенгенций, тест может работать некорректно\n",
      "chi2: 5.238095238095239\n",
      "p_value: 0.02209745528422724\n",
      "pearson: 0.5679618342470648\n",
      "chuprov: 0.6900655593423543\n",
      "kramer: 0.6900655593423543\n",
      "sks: 0.6900655593423543\n",
      "guttman: 0.0\n",
      "\n",
      "Между признаками в таблице сопряженности имеется зависимость на уровне значимости 0.05\n"
     ]
    }
   ],
   "source": [
    "alpha = 0.05\n",
    "bins = 2\n",
    "\n",
    "print(\"====================================\\n\")\n",
    "for test_name, test_func in zip(['pearson', 'spearman', 'kendall'], [pearson_test, spearman_test, kendall_test]):\n",
    "    stat, p_value, corr = test_func(df)\n",
    "    print(f'{test_name}_stat:', stat)\n",
    "    print(f'{test_name}_p_value:', p_value)\n",
    "    print(f'{test_name}_corr:', corr, '\\n')\n",
    "\n",
    "contingency_table = make_contingency_table(df, bins=bins)\n",
    "\n",
    "print(\"====================================\\n\")\n",
    "print(f'Таблица сопряженности для {bins} бинов:\\n')\n",
    "print(contingency_table, '\\n')\n",
    "\n",
    "chi2, p_value, chi2_norm = chi2_test(contingency_table, return_stats_norm=True)\n",
    "\n",
    "print('chi2:', chi2)\n",
    "print('p_value:', p_value)\n",
    "for chi2_norm_coef_name, chi2_norm_coef_val in chi2_norm.items():\n",
    "    print(f\"{chi2_norm_coef_name}: {chi2_norm_coef_val}\")\n",
    "print('guttman:', gutman_test(df))\n",
    "\n",
    "if p_value < alpha:\n",
    "    print(f\"\\nМежду признаками в таблице сопряженности имеется зависимость на уровне значимости {alpha}\")\n",
    "else:\n",
    "    print(f\"\\nМежду признаками в таблице сопряженности отсутствует зависимость на уровне значимости {alpha}\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
